from datetime import datetime

import pandas as pd
import numpy as np
import xarray as xr

import ttide
from ttide.t_tide import t_tide
from scipy.optimize import curve_fit


from ppgnss import gnss_time
from ppgnss import gnss_utils
from ppgnss import gnss_geodesy

from plot_baseline import read_seal_file, read_gamit_solution
import matplotlib.pyplot as plt
font_path = '/Library/Fonts/PingFang.ttc'  # 请替换为实际的字体文件路径

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
# 定义周期信号模型
def periodic_signal(t, A, phi):
    # 23小时55分55秒的周期对应的角频率
    omega = 2 * np.pi / (23 * 3600 + 55 * 60 + 55)
    return A * np.sin(omega * t + phi)

gamit_solution = "../baseline/solution.dat"

zihe_excel_file_path_02 = "../data/cczihe02.xlsx"
zihe_excel_file_path_03 = "../data/cczihe03.xlsx"
excel_file_path = zihe_excel_file_path_02

seal_down_file = "../data/sea_levels/downstream_7308.txt"
seal_up_file = "../data/sea_levels/upstream_3020.txt"

lat1, lon1 = 30.22863399, 120.73933035
lat2, lon2 = 30.22401975, 120.74423075

x1, y1, z1 = gnss_geodesy.blh2xyz(lat1, lon1, 0)
x2, y2, z2 = gnss_geodesy.blh2xyz(lat2, lon2, 0)

north, earth, up = gnss_geodesy.dxyz2neu([x2-x1, y2-y1, z2-z1], [x1, y1, z1])
alpha = np.arctan(north/earth)

df = pd.read_excel(excel_file_path)
xr_gg = read_gamit_solution(gamit_solution, 2022, 118)
# xr_seal_down = read_seal_file(seal_down_file)

# xr_seal_up = read_seal_file(seal_up_file)

# ['采集时间', '设备序列号', '东向坐标 (m)', '北向坐标 (m)', '天向坐标 (m)', 
# '东向坐标日变形速率 (mm/d)',  '北向坐标日变形速率 (mm/d)', '天向坐标日变形速率 (mm/d)',
# '平面日变形速率 (mm/d)', '三维日变形速率 (mm/d)', '差分期龄 (s)']
count = len(df["采集时间"])
time_from = np.datetime64(df["采集时间"][count-1])
time_to = np.datetime64(df["采集时间"][0])
dtime = np.datetime64(df["采集时间"][count-2])-np.datetime64(df["采集时间"][count-1])
nepoch = (time_to - time_from)//dtime+1

times = np.array([time_from + dtime*_ - 7.5*np.timedelta64(1, "h") for _ in range(nepoch)])
a_time_from, a_time_to = "2023-01-29 00:00:00", "2023-03-09 00:00:00" # 50 days
xr_gg = xr_gg.loc[a_time_from:a_time_to]

idx_from_s = np.argmin(df["采集时间"] >= a_time_from+ " 00:00:00")
idx_to_s = np.argmin(df["采集时间"] >= a_time_to+ " 00:00:00")
dus = df["天向坐标 (m)"][idx_to_s:idx_from_s]-np.nanmean(df["天向坐标 (m)"][idx_to_s:idx_from_s])
dns = df["北向坐标 (m)"][idx_to_s:idx_from_s]-np.nanmean(df["北向坐标 (m)"][idx_to_s:idx_from_s])
des = df["东向坐标 (m)"][idx_to_s:idx_from_s]-np.nanmean(df["东向坐标 (m)"][idx_to_s:idx_from_s])
err_max = 0.04 # 0.04
idx = np.logical_or(np.logical_or(np.abs(dus) > err_max, np.abs(dns) > err_max), np.abs(des)> err_max)
dus[idx] = np.nan
dns[idx] = np.nan
des[idx] = np.nan

dus = dus - np.nanmean(dus)
dns = dns - np.nanmean(dns)
des = des - np.nanmean(des)

dus_gg = xr_gg.loc[:, "u"].values - np.nanmean(xr_gg.loc[:, "u"].values)
dns_gg = xr_gg.loc[:, "n"].values - np.nanmean(xr_gg.loc[:, "n"].values)
des_gg = xr_gg.loc[:, "e"].values - np.nanmean(xr_gg.loc[:, "e"].values)

shape = (len(times), 5)
data = np.full(shape, np.nan)

xr_data = xr.DataArray(data, dims=["time", "data"], coords=[times, ['n', 'e', 'u', "v", "a"]])
xr_gg_data = xr.DataArray(np.full(xr_gg.values.shape, np.nan), dims=["time", "data"], coords=[xr_gg.coords["time"].values, ["u", "v", "a"]])
valid_times = np.array([np.datetime64(_) - 7.5*np.timedelta64(1, "h") for _ in df["采集时间"][idx_to_s:idx_from_s]]) 

v = des*np.cos(alpha) + dns*np.sin(alpha)
a = -des*np.sin(alpha) + dns*np.cos(alpha)

v_gg = des_gg*np.cos(alpha) + dns_gg*np.sin(alpha)
a_gg = -des_gg*np.sin(alpha) + dns_gg*np.cos(alpha)
# import sys
# sys.exit(0)

xr_data.loc[valid_times, "u"] = dus
xr_data.loc[valid_times, "n"] = dns
xr_data.loc[valid_times, "e"] = des
xr_data.loc[valid_times, "v"] = v
xr_data.loc[valid_times, "a"] = a


xr_gg_data.loc[:, "v"] = v_gg
xr_gg_data.loc[:, "a"] = a_gg
xr_gg_data.loc[:, "u"] = dus_gg


time_stick = ["2023-01-29", "2023-02-23", "2023-02-24", "2023-02-25", "2023-02-26", "2023-02-27"]
# xr_time1, xr_time2 = xr.align(xr_data.coords["time"], xr_seal_down.coords["time"])
# print(xr_time1)
time_resample = [np.datetime64(str(times[0])[:10]+" 00:00:00")+np.timedelta64(_, "h") for _ in range(500*24)]


xr_data = xr_data.loc[a_time_from:a_time_to]
dus = xr_data.loc[:, "u"].values
das = xr_data.loc[:, "a"].values
dvs = xr_data.loc[:, "v"].values

print("gg std: V: %10.6f  A: %10.6f  U: %10.6f"% (np.nanstd(v_gg), np.nanstd(a_gg), np.nanstd(dus_gg)))
print(np.nanmax(dus_gg), np.nanmin(dus_gg), np.nanargmax(dus_gg), np.nanargmin(dus_gg))
print(np.nanmax(v_gg), np.nanmin(v_gg), np.nanargmax(v_gg), np.nanargmin(v_gg))
print(np.nanmax(a_gg), np.nanmin(a_gg), np.nanargmax(a_gg), np.nanargmin(a_gg))

print("2h std: V: %10.6f  A: %10.6f  U: %10.6f"% (np.nanstd(v), np.nanstd(a), np.nanstd(dus)))
print(np.nanmax(dus), np.nanmin(dus), np.nanargmax(dus), np.nanargmin(dus))
print(np.nanmax(das), np.nanmin(das), np.nanargmax(das), np.nanargmin(das))
print(np.nanmax(dvs), np.nanmin(dvs), np.nanargmax(dvs), np.nanargmin(dvs))

enu_outfn = "enu_0129_0309.npo"
gnss_utils.saveobject(xr_data, enu_outfn)
width, height = gnss_utils.cm2inch(18), gnss_utils.cm2inch(10)
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(width, height), sharex=True, sharey=True)

plt.subplots_adjust(left=0.08, right=0.98, bottom=0.12, top=0.91, hspace=0)

coord_time = xr_data.coords["time"]
coord_time_gg = xr_gg_data.coords["time"]

xind = [(_-coord_time[0])/np.timedelta64(1, "D") for _ in coord_time]
xind_gg = [(_-coord_time_gg[0])/np.timedelta64(1, "D") for _ in coord_time_gg]

linewidth=1
markersize=3
h0, =axes[0].plot(xind, xr_data.loc[time_from:time_to, "u"], '-', linewidth=linewidth, label=u"2小时解")
h0d, = axes[0].plot(xind_gg, xr_gg_data.loc[:, "u"], 'o', c="r", label=u"日解", markersize=markersize)
axes[0].legend([h0, h0d], [u"2小时解", u"日解"], ncol=2, loc='upper left', bbox_to_anchor=(0.34, 1.35))

axes[0].set_ylabel(u"U (cm)")
h1, =axes[1].plot(xind, xr_data.loc[time_from:time_to, "v"], '-',linewidth=linewidth, label=u"2小时解")
axes[1].plot(xind_gg, xr_gg_data.loc[:, "v"], 'o', c="r", label=u"日解", markersize=markersize)

axes[1].set_ylabel(u"Y (cm)")
# xind = xr_delta.loc[time_from:time_to].coords["time"]
h2, =axes[2].plot(xind, xr_data.loc[time_from:time_to, "a"], "-",linewidth=linewidth, label=u"2小时解")
axes[2].plot(xind_gg, xr_gg_data.loc[:, "a"], 'o', c="r", label=u"日解", markersize=markersize)

axes[2].set_ylabel(u"X (cm)")
# plt.ylabel("Position Error (cm)")
y_lim_min, y_lim_max = -0.03, 0.03
axes[0].grid(axis='y', linestyle='-')
axes[0].set_ylim((y_lim_min, y_lim_max))

axes[1].grid(axis='y', linestyle='-')
axes[1].set_ylim((y_lim_min, y_lim_max))

axes[2].grid(axis='y', linestyle='-')
axes[2].set_ylim((y_lim_min, y_lim_max))

xticks = range(0, 50, 5)
xlbls = []
for idy in xticks:
    start_time = coord_time[0].values
    current_time = start_time + np.timedelta64(idy, "D")
    # print(current_time)
    lbl = str(current_time)[5:10]
    xlbls.append(lbl)
yticks = np.arange(y_lim_min, y_lim_max+0.001, 0.015)
ytklbls = [f"{_*100:.1f}" for _ in yticks]
plt.xticks(xticks, xlbls, fontsize=10)
plt.yticks(yticks, [None]+ytklbls[1:-1]+[None], fontsize=10)
plt.xlim(xind[0], xind[-1])
# plt.legend([h0, h1, h2], [u"闸外水位", u"闸内水位", u"内外位差"], loc='upper center', ncol=3, fontsize=10,bbox_to_anchor=(0.45, 1.26))
plt.xlabel(u"时间 (2023年, 月-日)", fontsize=10)
fig_filename = "figures/position_error.png"
plt.savefig(fig_filename, dpi=300)
plt.show()